A hierarchical phase space generator 
for QCD antenna structures 



Andre van Hameren|^ and Costas G. PapadopoulosQ 

Institute of Nuclear Physics, NCSR AijfioKgLTOc^, 15310 Athens, Greece 
g \ February 1, 2008 

O 

(N 



^3 



Abstract 



We present a "hierarchical" strategy for phase space generation in order 
. to efficiently map the antenna momentum structures, typically occurring in 

\ QCD amplitudes. 
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§ ■ 1 Introduction 
o 

g . The reliable description of multi-jet production at the LHC ^ will be an important issue. 
^ , This is not only related to the study of QCD in multi-parton final states but it is also very 
important in order to estimate several backgrounds for new physics effects. For instance, 
Qh! new unstable massive particles that decay to many partons may be discovered at the LHC 
r-| ! only when a reliable description of these final states is established. 

In this respect, apart from the problem of computing scattering matrix elements with 
^ ! many particles, also the efficient phase space generation is of great importance, because 
the scattering amplitudes in QCD exhibit strong peaking structures in phase space, which 
have to be taken into account by the generation algorithm. Flat phase space generators, 
like RAMBO [0, will not be adequate for this task. In the last years several methods to 
efficiently integrate the peaking structures of the scattering amplitudes have emerged, and 
have been used in several contexts 0. For instance, PHEGAS 0] is an example where 
an efficient, automated, mapping of all possible peaking structures of a given scattering 
process has been established. The algorithm is based on the "natural" mappings dictated 
by the Feynman graphs contributing to the given process, so that the number of kinematical 
channels used to generate the phase space is equal to the number of Feynman graphs. Using 
adaptive methods, like multi-channel optimization |^ and by throwing away channels that 
are negligible, we may end up with a few channel generator exhibiting high efficiency, as 
is indeed the case in n(+7)-fermion production in e^e~ collisions. In contrast, the QCD 
scattering amplitudes point towards the opposite direction: large number of Feynman 
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graphs which means large number of kinematical channels which, moreover, contribute 
equally to the result. 

A way out off this problem may be based on the long-standing remark that n + 2- 
gluon amplitude may be described by a very compact expression when special helicities 
are assigned to the gluons, which, combined with the leading color approximation, results 
to 

/ jST \n n+2 

5^|>[r = 8(^j (iV,2-l) J2 ^n+2(j>l,...,Pn+2) , (1) 

c ^ ^ l<«<j P(2,...,n+2) 

where Nc refers to the number of colors, 

An+2{Pl,--- ,Pn+2) ■= [ (Pl " P2) (^2 " Ps) " " " (Pn+1 " Pn+2) (Pn+2 " Pi) ] , (2) 

and the sum over all permutations of the 2"^^^ to the {n + 2)'^'^ argument of this function is 
taken, with the exception of those that are equivalent under reflection i n + A ~ i 0. 

SARGE 1^ is the first known example of a phase space generator that deals with the 
momentum structures entering the above expression, namely with (|^), known as antenna 
structures. The algorithm is based on the "democratic" strategy to generate the n body 
phase space, as is the case for RAMBO, and it makes use of the scale symmetry of the antenna 
to achieve the required goal. 

In this paper, we study the "hierarchical" strategy for phase space generation in order 
to efficiently map the momentum antenna structures. The idea is as follows. Using the 
standard two-body phase space (neglecting factors of 27r) 

d^2{P;si,S2;pi,P2) ■= d^PiS+{pl- si)d^P2 6+{pl- S2)S'^{P -P1-P2) , (3) 
we decompose the phase space as 

n n 

dMP: Pi • • • ,Pn) := ( n d% h{pl - a,)) 5^ ( 5^ p. - P) 

1=1 i=l 
= dSn-l Qn-l) 
X dSn-2 d^2{Qn~l'-, Sn-2'i Pn~l, Q n~2) 

X ds2d^2{Q3;(^3,S2;P3,Q2) d^2{Q2;cr2,ai]p2,pi) ■ (4) 

The task is to express the phase space in terms of the invariants Pi ■ pj appearing in the 
antenna structure (Q), so that, using a suitable mapping, we can construct a density that, 
apart from constant and soft terms, will be identical to this antenna structure. 

In the first section, we describe the basic building block of the algorithm, which is the 
expression of the two-body phase space in terms of the scaled invariants. In the second 
section, we demonstrate how this basic building block can be used in a sequential way 
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to produce the full antenna. Finally in the third section, some details concerning the 
numerical implementation of the algorithm as well as comparisons with known generators 
is given. The Appendices present all relevant generation algorithms from which the exact 
functioning of the generator can be reconstructed. 



2 The hierarchical antenna 
2.1 The basic building block 

To illustrate the idea we consider the generation of the 2-body phase space (0) when two 
massless antenna momenta qi.,q2 are given. The momentum P can be decomposed as 

p/^ = rgf + A'^ =: A'^ + 5f with r = P7(2P ■ gi) , 

so that the Sudakov parameterization of pi is given by 

= aiA^" + 6iPf + . 

where the variables ai and hi are given by 

Vi- Bi pi- Ai 

oi = w 7 oi 



Ai-Bi Ai-Bi 

The same can be done in terms of p2 and g2, and in the center-of-mass frame (CMF) of P, 
where P = {y/s , 0, 0, 0), cos(Z(g]^, gg)) = c and s = Vl — c^, one can choose 

= iv^(l, 0,0,-1) , v42 = |v^(l,0,-s,-c) , 
Bi = iv^(l, 0, 0, 1) , P2 = |v^(l, 0, s, c) , 

/Cl = (0,Xi,?/i,0) , /C2 = (0,X2,i/2C, -?/2s) . 

The phase space can now be completely expressed in terms of ai and 02, leading to 
with 

n(ai, 02) = 4s^[(l - a2 + S2 - Si)a2 - 82] 

- [(1 - 2ai - si + S2) + (1 - 2a2 - h + 82)0]^ , (6) 

where si^2 = 81^2/ s, and where B is the step function. In terms of Lorentz invariants, the 
parameters oi and 02 are given by 

qi-Pi P2-q2 
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So in order to obtain a two-body phase space with a density which depends on the invariants 
Oi, a2 following some given function /(cti, 02), one has to generate Oi, 02 following a density 
proportional to /(oi, 02) x n(~^/^)(ai, 02) in the region where n(ai, 02) > 0, and construct 
the momenta following the Sudakov parameterization. Explicitly, the direct construction 
is given by 

P?^(s + Si-S2)/(2v^) , 

pI^ Pi~ Vsai , 

pI^ iiVs ~Pi- Vsa2) + cpl )/s , 

pi^eiip',r-s,-ipir-ipiry/^ , 

where e should be a fair random variable which can take values +1 and —1. For more 
details about this procedure, we refer the reader to the Appendix 



2.2 Antenna generation 

In the hierarchical/sequential approach, the generation strategy proceeds through a se- 
quence of two-body phase space generations following the decomposition (||). At each 
two-body generation, one final-state momentum pk is generated, together with the sum 
Qk~i of the remaining final-state momenta to be generated. This suggest to label the 
momenta in a way opposite to the order of generation, so first p„, Qn-i are generated, then 
Pn-i,Qn-2 and so on. The starting point is the CMF of the initial momenta qi and q2 
with Qn = qi + q2 being the overall momentum. The CMF of momentum Qk we denote by 

(k) (k) 

CMFfc. The pair pk,Qk-i is generated by generating variables a\ ,03 and constructing 
the momenta as described before. These variables are now equal to 

(k) Pk+i ■ Pk 1 (k) q2 ■ Qk-i 

' — and a\ ' — 



Pk+i - Qk q2- Q 



k 



This happens in CMFfc, so in order to obtain Pk,Qk-i, the constructed momenta have to 
be boosted such that (a/Q^, 0, 0, 0) is transformed to Q 



k- 

We would like to generate the momenta following a density that is proportional to 

An+2{qi,Pn--- ,Puq2) = [{qi-Pn){Pn-Pn-l) ■■■{Pl-q2)]~^ ■ (8) 

Since the integrand is infra-red singular, a cutoff on the invariants is necessary. Therefore, 
we define a symmetric matrix cjjj which encodes the restrictions on the momenta through 

k 

(^ii = (^i= p1 , y/cTiCTj < cTij < Pi ■ Pj and Ek := ^ o"* • (9) 



i=l 



Before we proceed, we do three observations. Firstly, we have 

Pk+i ■ Qk = {Ql+i -Ql- Pfc+i)/2 . 
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Secondly, we have 



— log =: gk+i{sk) , (10) 



{sk+i — — Sk) {sk — Sk) dsk \Sk+i — o'k+i — Sk 
and thirdly, we can write 

1 f ^ 1 \ 1 

A„+2(gi,P„. . . ,Pl,g2) = TZ^iSn - Sn)iqi ■ Qn){q2 ' Qn) TT ^fc(gfc-l) ~Ti^ ' 

with Sn = Qn: Pn+1 = and Qi = pi- These observations suggest that the phase space 
generation 

dSn-2 gn-l{Sn-2) rfa^""^^ ^^4""^^ nj^^n^ 9 (n(„_i) ) 



a 



ds2g?,{s2) daf^-^ ^4^^ nj3)^/^^e(n(3)) 

daf^^ daf^^ U^^/'^mm) ' (11) 
a\ ay 

will lead to a density for the momenta that is proportional to An+2- Three variables are 
generated in each CMFfc, namely Sk-i, af^^ and a''2^. Just as the integration of Sk-i (|To|), 
also the integration of a^i \ a'^'^ results in a volume factor that depends on the corresponding 
variables generated in CMFfc+i. As we will show in Appendix 0, however, these factors 
are logarithmic functions of their arguments and exhibit a non-singular behavior, and we 
call them soft factors. The total actual density will therefore be the product of n — 1 soft 
factors times the antenna structure under consideration. 

In the end, we want to generate all permutations in the momenta of (^. Those for 
which qi and q2 each appear in two factors (none of which is qi- q2) cannot be obtained by 
simple re-labeling. In order to obtain these, we observe that they can be decomposed into 
two antennas, namely 

(91 5 Pm, Pm-l, ■ ■ ■ , P2, Pi, (I2) X An~m+2{(l2, Pn, Pn-l, ■ ■ ■ , Pm+2,Pm+l, Ql ) (12) 

and each of these can be generated after the decomposition, 

d^n , Pi ■ ■ ■ , Pn ) ds^ dSn—m d^2 i^Qm ^m, ^n—m, Qm, Qn—rri) 
X d^rn{Qni;Pl, ■ ■ ■ , Pm) d^ 

n—m {Q n~m ■,Pm+l,--- ,Pn) ■ (13) 

In order to combine the two sub-antennas to the required antenna structure, we have to 
take into account in the first decomposition a density that is proportional to 

©( y/^n 'sf^m \/ Sn—m ) 



{qi ■ Qm){qi ■ Q n—m ){q2 ■ Qm){q2 ■ Q n—m) ^m ^n- 



(14) 
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The case of m = 1 is still special. Then, the first step in ([TT|) should be replaced by 

and the rest of the sequence should go on with the replacement of p„ by Pn+i at the second 
step. For the density g' we refer to Appendix |^. There, we have collected all integrals and 
generation algorithms related to the antenna generation described so far. 

2.3 Open antennas 

As it will be clear from the numerical analysis presented in the next section, the soft 
factors appearing in the description of the QCD antenna contribute to a certain extent to 
the variance of the Monte Carlo integration. There is an alternative approach, that still 
follow the hierarchical/sequential generation strategy, and give better results. It based 
on the observation that the production of an 'open' antenna structure, namely one where 
the last product q2 ■ pi is missing, is simpler, since it can be constructed without using 
the variables a'^\ They can be replaced by flatly generated azimuthal angles, to that 
da^^U^^'^l'^'^^ Q{Il(^k-j) dif^'^K The basic decomposition therefore becomes 

dSn-ign{Sn-l) dif^'^'^ da^^^ X dSn-2gn-l{Sn~2) dif^'^''^^ da!{~^^ X • • • 

• • ■ X ds2g3{s2) dip^^^ daf^ X daf^^ . (16) 

This way, we will get the antenna density (H) without the factor pi ■ q2 in the denominator, 
which is the reason why we call this an 'open' antenna. A 'closed' antenna can be obtained 
using the fact that, by combining two open antennas, one can choose for another factor from 
the antenna string to be missing. Then, a multi-channeling procedure can be performed 
with these different choices, leading to a density that is, roughly speaking, proportional to 

{qi ■ Pn) + {Pn ■ Pn-l) H ^(^3-^2) + (P2 " Pi) + jPl " ^2) 

{qi ■ Pn){Pn ■ Pn-l) " " " (PS " P2){P2 " Pl){pi " ^2) 

To get the different choices, a first splitting of Qn into Qm and Qn-m has to be performed, 
after which open antennas are generated from each of these, one with qi and the other 
with q2 as initial momentum. For details, we refer to Appendix 0. 

3 Results 

In this section, we present results obtained by SARGE and HAAGQ, the program that imple- 
ments the hierarchical algorithm described before. In order to be as general as possible, 

■^HAAG stands for: Hierarchical AntennA Generation. 
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the only cut we apply is 



{Pi + Pjf > So , 

where 7^ j) runs from 1 to n + 2 where n is the number of final-state particles. 

Unless explicitly mentioned differently, we use sq = 900 GeV^ and the total energy ys = 
1000 GeV. Moreover, all particles are assumed to be massless in order to compare with 
SARGE, with which only massless particles can be treated. 

As it was mentioned in the introduction, we are interested in integrating sums of QCD 
antenna structures (^. We start by considering the simplest case, namely integrating the 
function 

{Pl ■ P3){P3 ■ Pi){Pi ■ P2){P2 ■ Ps) • • • {Pn+2 ' Pi) (17) 

that corresponds to a given permutation of the momenta, namely (1, 3, 4, 2, 5, . . . ,n + 2). 
In Table |I] we give the results for SARGE, HAAG with open antenna generation, and HAAG(C) 
with closed antenna generation. In all three codes the same single channel, corresponding 
to (p!?!) , has been used in the generation. A^gen and A'acc are the number of generated and 
accepted events, and by / we define 

f ■= — 

J ■ J2 ' 

where V2 is the quadratic variance and / is the estimated integral. / is clearly a measure 
of the efficiency of the generator. Moreover e, defined as 

<w > 
e := , 

^max 

is the usual generation efficiency related for instance to 'unweighted' events in a realistic 
simulation. The results agree well, and exhibit the fact that the generated densities of the 



jets 


algorithm 


-Agen 


A^acc 


/ 


A/ 


/ 


e{%) 


4 


SARGE 
HAAG 
HAAG(C) 


1 X 10^ 
5 X 10^ 
5 X 10^ 


34853 
31193 
28366 


.251 X 10'9 
.260 X 10-9 
.256 X 10-9 


.734 X 10-11 
.280 X 10-11 
.252 X 10-11 


85.9 
5.75 
4.84 


0.34 
1.77 
4.22 


5 


SARGE 
HAAG 
HAAG(C) 


2.5 X 10^ 
6.5 X 10^ 
6.5 X 10^ 


30960 
29855 
24345 


.438 X 10-1° 
.442 X 10-1° 
.441 X 10-1° 


.153 X 10-11 
.640 X 10-12 
.706 X 10-12 


307 
13.6 
16.7 


0.23 
1.02 
1.04 


6 


SARGE 
HAAG 
HAAG(C) 


1 X lO'' 
1.2 X 10^ 
1.2 X 10^ 


28383 
32070 
25040 


.487 X 10-11 
.487 X 10-11 
.485 X 10-11 


.164 X 10-12 
.658 X 10-13 
.886 X 10-13 


1141 
21.9 
40.1 


0.21 
1.48 
0.69 



Table 1: Results for the single-channel integration/generation. 



generators the hierarchical type are much closer to the integrand. Moreover, the closed 
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antenna algorithm HAAG(C) becomes less efficient compared to the open one as the number 
of particles increases. The same picture is reproduced for an arbitrary permutation. 

For a realistic QCD calculation, the integrated function may be approximated by a sum 
over permutations. Therefore, an efficient generator has to include all possible channels, 
where each channel corresponds to a given permutation of the momenta. In that case, a 
multi-channeling optimization procedure can be applied, which is incorporated in HAAG. 
In order to study the efficiency of this optimization we consider the same integration as 
before, but with all channels contributing to the generation and allowing for optimization. 
In this optimization procedure, we discard channels that contribute less than a certain 
pre-determined fraction to the set of available channels. It is expected, of course, that in 
end the right permutation will be 'chosen' by the optimization. This is indeed the case and 
the results are presented in Table ^ We see that the optimization results to a picture close 
to the one obtained with the single channel generation, with some noticeable improvement 
in the case of SARGE. We also include results with SARGE.n, a slightly different version, 
described in Appendix 0. 



jets 


algorithm 


Ngcn 


A^acc 


/ 


AI 


/ 


e{%) 




SARGE 


1 X 10^ 


52516 


.262 X 10" 


-9 


.294 X 10" 


-11 


12.6 


1.29 


4 


SARGE . n 


1 X 10^ 


46529 


.260 X 10- 


-9 


.298 X 10- 


-11 


13.2 


1.55 




HAAG 


5 X 10^ 


34293 


.257 X 10- 


-9 


.210 X 10- 


-11 


3.36 


4.28 




HAAG(C) 


5 X 10^ 


29736 


.259 X 10" 


-9 


.227 X 10" 


-11 


3.84 


3.91 




SARGE 


2.5 X 10^ 


32315 


.422 X 10" 


-10 


.106 X 10" 


-11 


159 


0.44 


5 


SARGE . n 


2. X 10^ 


30994 


.440 X 10" 


-10 


.807 X 10- 


-12 


67.2 


0.83 




HAAG 


6.5 X 10^ 


31063 


.444 X 10- 


-10 


.503 X 10- 


-12 


8.32 


1.17 




HAAG(C) 


6.5 X 10^ 


24179 


.436 X 10- 


-10 


.593 X 10- 


-12 


12.03 


1.84 




SARGE 


1 X 10^ 


29138 


.476 X 10" 


-11 


.145 X 10" 


-12 


933 


0.45 


6 


SARGE . n 


1 X 10*^ 


35445 


.492 X 10- 


-11 


.109 X 10- 


-12 


492 


0.25 




HAAG 


1.2 X 10^ 


33278 


.483 X 10- 


-11 


.595 X 10- 


-13 


18.2 


1.19 




HAAG(C) 


1.2 X 10^ 


24126 


.471 X 10" 


-11 


.749 X 10" 


-13 


30.3 


1.21 



Table 2: Results for the all-channel generation with optimization. 



As is the case for any multi-channel generator, a computational complexity problem 
arises when the number of channels increases. For instance, in our case we are facing a 
number of |(n -f- 1)! channels! On the other hand, it is also clear that the channels we are 
considering have a large overlap in most of the available phase space. It is therefore worth 
to investigate the dependence of the integration efficiency on the number of channels used. 
This is presented in Table ^ where the full antenna 

[ (Pl ■P3)(P3 ■P4)(P4 ■P2)(P2 -Ps) • • • (Pn-f2 -Pl) (18) 

P(2,..,n+2) 

is integrated, using a number of channels that has been selected on a random basis. We see 
the rather interesting phenomenon that a decent description can be achieved with a much 
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smaller number of channels. Variations of this technique of using only subsets of channels, 



# channels 


2520 


1500 


1000 


500 


200 


50 


10 


/ 


5.33 


5.37 


5.48 


5.72 


6.14 


11.6 


84.7 


'acc 


26630 


26521 


26437 


26676 


27009 


27190 


27205 


e{%) 


11.2 


13.1 


11.6 


7.1 


7.5 


1.7 


0.28 



Table 3: All-channel integration with subsets of channels for generation. 



for example choosing another subset after each step of multi-channel optimization, lead to 
the same picture. 

The complete results of the integration of the full antenna are presented in Table ^. 
We see that HAAG has a much better / factor than SARGE. On the other hand the e exhibits 



jets 


algorithm 


-^gcn 




/ 


AI 


e{%) 


/ 


4 


SARGE 
HAAG 


1 X 10^ 
6 X 10^ 


47483 
42019 


.166 X 10-^ 
.167 X 10-^ 


.115 X 10-9 
.810 X 10-1° 


4.21 
12.01 


4.8 
1.4 


5 


SARGE 
HAAG 


3 X 10^ 
1.2 X 10^ 


39095 
55234 


.176 X 10-^ 
.177 X 10-^ 


.162 X 10-9 
.856 X 10-1° 


3.27 
7.53 


25.6 
2.7 


6 


SARGE 
HAAG 


1.5 X 10^ 
1.8 X 10^ 


44529 
47911 


.157 X 10-^ 
.161 X 10-^ 


.135 X 10-9 
.905 X 10-1° 


2.95 
7.15 


109 
5.7 


7 


SARGE 
HAAG 


1 X 10^ 
3.6 X 10^ 


47766 
45599 


.123 X 10-^ 
.123 X 10-^ 


.988 X 10-1° 
.241 X 10-1° 


3.02 
5.11 


642 
13 


8 


SARGE 
HAAG 


1 X 10^ 
1 X 10^ 


53560 
49206 


.784 X 10"^ 
.789 X 10-8 


.554 X 10-1° 
.496 X 10-1° 


3.29 
6.30 


4998 
39 



Table 4: Results for the all-channel integration. 



a less dramatic effect. This is related to the fact that SARGE generates a phase space that 
is much larger than the one defined by the cut on sq. In that sense, if the main time 
consumption in a given computation is spent over the evaluation of the integrand (matrix 
element squared), it is more fair to compare the square of the estimated expected error, 
normalized by the number of accepted events A^^acc- In that case we see that HAAG is still 
2-3 times more efficient, and if we consider a smaller cut, namely ^/sq = 10 GeV, this gain 
goes up to an order of magnitude (Table 

For a realistic calculation of the cross section of a QCD process, one may assume that 
the time it takes to perform one evaluation of the integrand is much larger than the time 
it takes to generate one accepted event and to calculate the weight. This means that the 
computing time is completely determined by the number of accepted events A^acc- We 
introduce 

iVacc/ 

N 

^ 'gen 
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jets 


algorithm 


N 

^ ' gen 


N 

^ 'acc 


/ 


A/ 


e{%) 


/ 


4 


SARGE 
HAAG 


1 X 10^ 
6 X 10^ 


60986 
46763 


.364 X 10-^ 
.366 X 10-6 


.548 X lO-s 
.235 X 10-s 


0.631 
4.34 


22.7 
2.47 


5 


SARGE 
HAAG 


2 X 10^ 
1 X 10^ 


43150 
56034 


.619 X 10-6 
.643 X 10-6 


.165 X 10-^ 
.465 X 10-* 


0.29 
1.84 


142 
5.23 


6 


SARGE 
HAAG 


1 X 10^ 
1.4 X 10^ 


67811 
51983 


.114 X 10"^ 
.111 X 10-5 


.257 X 10-^ 
.883 X 10-8 


0.28 
2.50 


502 
8.83 


7 


SARGE 
HAAG 


5 X 10*^ 
2 X 10^ 


84391 
44015 


.186 X 10-5 
.192 X 10-5 


.346 X 10-^ 
.177 X 10-^ 


0.176 
2.24 


1723 
16 


8 


SARGE 
HAAG 


5 X 10^ 
5 X 10^ 


175541 
58874 


.354 X 10-5 
.350 X 10-5 


.517 X 10-^ 
.289 X 10-^ 


.119 
1.65 


10618 
34 



Table 5: Results for the all-channel integrationwith sq = 100 GeV^. 



as a measure of the computing time. For a realistic calculation, one has to multiply this 
number by the evaluation time of the integrand, and devide by the square of the relative 
error one wants to reach. Figure 1 shows this quantity as function of the number of 
produced partons using the data of Table ^. According to this graph, a calculation with 
SARGE would take 10 times longer than the calculation with HAAG. 



40 
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Figure 1: N^ccf /Ngen (a measure of computing time) as function of the number of produced 
partons. 



Finally, Figure 2 shows the dependence of the result on the value of the infrared cut- 
off y/sQ and the number of produced partons. The function clearly exhibits a negative 
power behavior for sq. Moreover, the curve becomes steeper with increasing number of 
jets, suggesting that at least the leading power is related to the final-state multiplicity. 
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Figure 2: The integral of the full antenna as function of ^^/so for different values of the 
number of produced partons. 

4 Conclusions 

HAAG exhibits an improved efficiency compared to SARGE for multi-parton calculations. 
It is also more powerful in describing densities where a partial symmetrization over the 
permutation space is considered. Finally, HAAG makes no fundamental distinction among 
massless and massive particles, so it can be used for an arbitrary multi-partonic process. 
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Appendices 

The following appendices contain details about the generation of the various random vari- 
ables necessary to build phase space with the desired density. The techniques used to 
achieve this are inversion, rejection and multi- channeling. For details about these tech- 
niques, we refer to 0] . We only want to mention that inversion is applicable if one has an 
analytic expression (with reasonable complexity) of the inverse of the indefinite integral 
of the density. Rejection can be applied if one knows a function which is strictly larger 
than the density and is proportional to a density one is able to generate. The efficiency 
is given by one divided by the integral of that function. Multi-channeling can be used if 
the density can be written as the weighted sum of densities, each of which one is able to 
generate, and if the weights are positive. 



A Closed antenna 

In this appendix we present all relevant algorithms for the generation of the closed antenna. 
More specifically, in the following we describe the generation of the three variables Sk-i, 

(k) (k) 

a\ and 03 that are needed to describe the antenna at the k—th CMF. Most of the time 
the superscript (k) is omitted for convenience. 
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Al Generation of Sfc-i 

Since the indefinite integral is given in (|TUp, we see that Sk~i can be generated with inver- 
sion. The hmits are given by 

k k 

Sk - Ak-i > Sk-i > Afc-i with Afc = IJk + ^ (Tij , Ak = cr^+i + 2 ^ (Tk+i,i , (19) 

and the weight factor is 

, I Sk ~ Sk-i — Ak-i\ , I Sk — (Jk — Afc_i 
log X + log 



Ak-1 — Ok J V Ajt_i — Ek 

In the case of m = 1 (|15D we simply have 

dsk-i ^-^ — = log(sfc_i - Sk-i) , (20) 

Sk-l — ^k-l 

SO that the weight factor is 

log {sk - Afe_i - Ek-i) - log (Afe_i - Ek-i) . 

A2 Generation of ai, a2 

The general integral corresponding to the generation of ai, 02 is given by 

e(n(ai,a2))e(ai-af^)e(a2-af ) 

aaiaa2 , (21) 

aia2A/n(ai, 02) 

where 11 is defined in (^, and possibly necessary infra-red cut-offs. We shall 

analyze this integral by first integrating over the a2-variable, and then over the oi-variable. 
The generation has then to be performed in the opposite order: first Oi and then 02. As 
we shall see, the inclusion of the cut-off on ai does not lead to complications, but the 
inclusion of both cut-offs does. For that case, we see two solutions. Firstly, we can replace 
the integral by 

, , e(n(ai,a2))e(ai-af^) 
aaiaa2 



ai{a2 + ai, 02 j 

where h is related to This, of course, changes the actual density with which ai and 
a2 are generated, but may, for small h, still be considered suitable for the desired antenna 
structure. Also, the final antenna generator will cover phase space less efficiently (will 
generate 'to much' phase space), simply because less cuts are included analytically. In the 
second solution, we write the integrand as the sum 

e(n(ai,a2))9(ai -gf^) _^ 9( n(ai, 02) ) e( as - 4°^ ) 

0-1, 0-2 

) (ai + 02)02^/11(01,02) 
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which is exactly equal to the original integrand on the phase space for which both cuts 
are included. Both integrands can be integrated analytically, so that the multi-channeling 
procedure can be applied to generate their sum. For this solution, the only problem is 
that, again, 'to much' phase space is generated. 

A2.1 The a2-VELriable In order to integrate over the 02- variable, n(ai,a2) is more 
conveniently written as 

n(ai,a2) = 4(a^(ai) - a2)(a2 - 02 (ai) ) , (22) 

with 

«2 — o H h 1 - 2ai c 

2 \ Sk \ Sk / / 



±JaJl-a,-'-^^^)-^Y' . (23) 



Sk / Sjfc / 

Then, the general a2-integral is given by 



/ 

J a. 



'2 

da2 

2 



(a2 + /i) V {at - 02) (02 - 0^2 ) 



with /i = if no cut-off on 02 is desirable, and h related to af^ or h = ai, depending on 
the solution mentioned above if a cut-off is desirable. By substituting a2 <— a2 — h, this 
integral can be written as 



1 -2 



/ = = 

Jd^ 02 v (^2 ~ ct'2){a2 — 0,2) ya^a^ 



arctan 



0-2 - Q2j 

at{a2 - 0,2) 



ai a2 

with af — af + h. The explicit indefinite integral shows that the variable 02 can be 
generated by inversion. 

A2.2 The ai- variable We start this section by mentioning that, in the case = 
(Tk+i 7^ 0, the variables 

(k) Pk ■ Pk+l 
a] ' 



^ Qk ■ Pk+l 

can be expressed in terms of the 'massless' af^^ defined in terms of the 'long' component 
pi+i of Pk+l in the CMFfc: if = (pO+^, Pk%+i) and p[^}^ = {pl^^, -Pk%+i) with 

^ |Pfc+il ^/X{sk+l,ak+l,Sk) 

Pk ■= — — = , l^4j 

Pk+l ^k+l — — Sk 
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then = + ^Pi+i ' and 

Since the relation between and a^'^^ is hnear, the two-body phase space is still expressible 
in terms of a^'^^ scaled by jSk- The generic oi-integral is given by 

:iax ^ 

dai , (25) 

(ai + /ii) V ( ('^i) + ^(oi) )( 02 (fli) + /i(ai) ) 

with as defined in (P^D, his a. constant (for this integral) or h[ai) = ai. The kinematical 
integration limits, coming from the requirement that (ai) are real, are given by 



2 Y Sk ]l \ Sk Sk 

In the massless case we get = 0,ai = 1 — Sk-i/sk and we have to impose a lower bound 
on ai given by a^°^ = (Jk+i,k/Pk+i ■ Qk- In all cases for the form of h{ai), the ai-integral is 
of the type 

1 



dai- — - with f (ai) := \/ + 2vai + w"^ , v"^ < , (26) 

(ai + /ii)/(ai) 

and the indefinite integral is given by 

1 1 f ai + hi + f{ai) - f{-hi) 

dai 7 — - = — — — log 



(ai + /ii)/(ai) /(-/ii) " V«i + ^1 + /(«i) + /(-^i), 
which is analytically invertible. The definite integral (times tt/A from the a2-integral and 
a factor l/Pk if Pk+i is massive) gives the weight factors in the generation of ai, a2 for the 
closed antenna. 



A2.3 The ai-variable in the case m = 1 This refers to (ITH), where we have c = — 1 
and a2 = ai + fi with 

^ := ^ 1 , so that / dai —jz r-, ■ — rrz r (27) 

Sk Amin ai{l-ai){ai + fi-ai) 

is the integral to be performed. Since the integrand is equal to 

111 1 
+ r + 7 rrz r + 



ai{ai + jj,) ai(l— /i — fli) (oi + /i)(l — Oi) (1 — ai)(l — /i — ai) 

we see that the generation of ai can be done easily using the multi-channel technique with 
four channels with weight 

= -l2) 0)7 (2) , - log (2) , 

where 
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i 


di 


ei 




(2) 


1 





fi 


+ 


+ 


2 





1 — /i 


+ 




3 


/i 


1 


+ 




4 


1 


1 — ;U 







The soft factor will simply be equal to Ylt=i "^^^ integration over the azimuthal angle, 
replacing the a2-integration, gives an extra factor of 27r. 



A3 Antenna split 



With qi oc (1, 0, 01) and g2 oc (1, 0, 0, 
the form 



-1), the two-body phase space integral (|14|) assumes 



dsids2 dQz 
with the energies 

Ei{si, S2) ■■ = 

M 



(0) 



S1S2 ( E,isu S2)2 - )( E2(si, - Ql 



f28) 



S + Si - S2 



-^2(51, S2) := y/s - Ei{si, S2) 



and where s\ 2 are the sums of the matrix elements aij (^) corresponding with the mo- 
menta in two antennas to be generated. Qz is integrated between the kinematical limits 
E\ — and can be treated in a way similar to the one described in the previous 
paragraph, by mult i- channeling over four channels with 



i 


di 


di 




(2) 

91 


1 


El 


E2 


+ 


+ 


2 


El 


E2 


+ 




3 


El 


E2 


+ 




4 


El 


E2 







The final weight is the sum of the channel-weights, divided by the (dimensionful) factor 
4E1E2. The generation of Si, S2 is a specific case of a more general problem described in 
Appendix [B2^. 



B Open antenna 

The three variables needed to describe the open antenna are Sk-i, oli^ and Lp^^\ In each 
CMFfc, the (^'^'^^-variable should be generated with uniform distribution between and 2-7?, 

ik) (k) 

and the a\ -variable should be distributed following l/cti between 

Sfc + (Tfe — Sk-l , Cfc, Sk-l) 



max 



Sk + CTfc 



2sfc 

Sk~l 



Pk^/Xi 



S, CTfc, Sk-l, 



2s, 



0'k+l,k 
Pk+1 ■ Qk 
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with I3k as defined in (24). The normahzed density for these generations is equal to 

with Lk = log{af^]_/a^i^_) . 



2/?fc _ A/A(sfc+i, CTfc+i, Sk) ^ .^^ ^ _ (fc) , (A;) 



TTLk aP T^Lk {pk+i ■ Pk) 
This suggests to use 



9k+i{sk) ■-- 



(sk — Sk) a/ A(sfc+i, o"fc+i, Sfc) 

■Sfc — ^fc + a/ (Tfe+i, Sfc) — \/ X{Sk+l, (Jk+l, Sk 




Sk — Sk + a/ X{Sk+l, (Tfe+l, Sjfc) + A/A(Sfe+i, (Tfc+i, Z'fc 

instead of ( p!OD for the case that 7^ 1. The logarithm contributes to the weight as a soft 
factor again. For small values of the squared masses ak, the factor in the numerator will 
be canceled by Sk+i — Sk+i in the denominator of gk+2{sk+i)- Just as in the case of the 
closed antenna, we end up with one remaining, and desirable, factor S2 — S2 = 2(p2 ■ Pi) in 
the denominator of the open antenna density, which cannot be achieved by the generation 
of the a['^''-variables. 

Let us denote the soft factor coming from the Sfc-generation by 



Gk+i '■= / dsk gk+iisk) , 



then the complete density resulting from the decomposition ([Tq) is 

dDo{Qn]Pn, ■ ■ ■ ,Pl) _ {qi ■ Qn) a/ A(s„, (T„, S^-l) 



d^n{,Qn]Pn, ■■■ , Pi) (QI ' Pn){Pn ' Pn-l) " " " (PS ■P2)(P2 " Pi) 

with soft factor 



Aak-iSk-2 

{A7t)^'^GsLsL2 y GkLk V ^ {sk-i - Sk-,) 

k=n " 



The factor pi ■ q2 is missing in the denominator of the density, which is the reason why 
we call this an open antenna. Other open antennas can be obtained by starting with a 
decomposition of Qn into two momenta, from each of which open antennas of the above 
type are generated, with initial momentum for the one, and q2 for the other. 

In order to digress about this procedure, let us extend the labeling a bit. With a set 
{/i, I2, ■ ■ ■ , In} of n non-equal labels, we write 

k 



m=l 

and so on. Now take n = M + N and let 



{h,... ,In}:={M + 1,M + 2,... ,M + N} , {J,, Jm} := {M, M - 1, ... ,2,1} . 
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We introduce pif^_^_i — qi, PJm+i — Q2, and the decomposition 

ri) Slj^ , Sjf^ , Qlj^ , Qjm ) 

X dDo{Qi^;pj^,pi^_^, . . . ,pi^) dDo{Qj^;pj^,pj^_^, . . . ,pj^) , 
which produces the density 

dDM{Qn;Pn, ■■■ ,Pl) _ 2(5^ Bj^jBij^ {pm+1 ■ Pm) 

X 



d^„{Qn]Pn, ■ ■ ■ ,Pl) Tr^/X{Ql,Sj^,Si^) (qi ■ Pn){Pn ■ Pn-l) ■ ■ ■ {P2 ■ Pl)(Pl ■ q2) 

X {qi ■ Qi^){q2 ■ QJM)^J Ksin,(^In, ^In-i)KsJm,(^Jm, ^Jm-i) ■ (29) 

In order to cancel the 'undesirable' factors on the second hne, we need to take care of the 
generation of s/^, sj^^, Qi^, Qj^. 

Bl Generation oi Q , Q 

Since qi oc (1, 0, 0, 1) and 52 oc (1, 0, 0, —1), we can write 

I1 ■ Qin 1 QJm ■ Q2 , ., 1 Sjj^ - Si^ 

oi = — — and — — ai + II with /i — , 

qi ■ Qn Qn ■ ?2 Sn 

SO that the generation of an azimuthal angle (/? between and 27r with the uniform distri- 
bution, and the generation of ai with a density proportional to 

^ between ^ + sj, - sj^ ± y^X{s^,s,,,sjJ 



aiiai + fji) 2sn 
leads to the total density 

dUQn. Q.., . ^^^^^^^^ . . 

B2 Generation of S/^^ , Sj^ 

B2.1 If Af = 1, then sj^^ — (Tj^i and we only need to generate s/^ with a density 
proportional to 



d 



log (^si^ - ai^ - Ei^_^ + ^\{si^,ai^,Si^_S) 



This density cancels the corresponding factor in the total antenna density. Something 
similar can be done if A?" = 1. 
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B2.2 If A?" > 1 and M > 1, then both s/^ and sjj^^ have to be generated, in the region 
where ^/s^— y/sj^ — ^/sj^j > 0. This is far more comphcated than the previous case, and we 
restrict ourselves to a density with a denominator proportional to {sj^ — Ej^){sjj^j —Ej^j). 
Because ct7^Z'/^_j and aj^jEjj^j_^ may be considered small, it still cancels the factor 



in (|29|) . We shall write Si, S2 instead of s/^, sj^^ from now on, and denote m := ^/s^ , 

N M 

ml := Ui^ , cl:= aj^j^ and ml := Tj,, , := ^ (jj^ . 

X,L=1 K,L=l 

We choose first to generate si, and then S2, so that the integral, corresponding with the 
generation, becomes 



2 Si - mj J ^2 S2-ml 

The S2-integral is simple, and shows that S2 can easily be obtained by inversion. After 
integration over S2, the Si-integral becomes 

/ — - — 2 [log (("^ ~ v^)^ ~ ^l) - log {4 - ml)] . 

The si-variable distributed following this integrand can be obtained with high efficiency by 
rejection from the density proportional to l/(si — ). The total integral can be calculated, 
and is given by 



p,e=±l 



where 



Li ( "^ + "^^ 1 + log ( ^ + P^^-^A log ( ^^+"^^ \\"^ 
m + pm2 + emi J \ m J \m + pm2 + emi J J ^^^^^ 

-iog(±:p) log f -f^^' 7 



Li(x) := J2 



n=l 



C SARGE 

In this Appendix, we give a short review about SARGE and present the adaptations applied 
in SARGE. n. We need to start with the establishment of some notation. 



19 



Tip is a Lorentz transformation that boosts momentum p to (a/j?, 0, 0, 0). TZp is a 
Lorentz transformation that rotates momentum p to 0, 0, The standard repre- 
sentation of a unit vector in terms of parameters z G [—1,1] and (p G [0, 27r] is denoted 
by 



n[z. 



(p) := (Vl — -2^ Sincyj, Vl — -2^ COS(y9, 



Pn C [— 1, 1]" is the n-dimensional polytope, which consists of the support of the indicator 
function 

n n 

'dp„{Xi,X2, . . . ,Xn) ■■= ]^0(1 - \Xi\) Y\_ - - ^k\) ■ 

i=l j,k=l 

The algorithm for the generation of n final-state momenta with a single antenna structure 
without initial-state momenta as given in is 

Algorithm C.l 

1. generate two massless momenta qi,qn, back-to-back; 

2. generate {xi, . . . , X2n-4) £ P2n-4 and {(p2, • • • , Pn-i) G [0, 27r]"^^, all uniformly dis- 
tributed; 

3. for i = 2, . . . ,n — 2 construct qi following (with xq := 0) 

h ^ qi-i + q-n , ^2 ^ 'Hb^qi-i , 



V ^ \l 7, (6 + 6) , , V ^qU^ n{z, pi) , 

« + u 



4. put Pi ^ uHqqi for i = 1, . . . , n, where Q = X]r=i ^^"^ ~ \f^jQ^ ■ 
The density with which the generated momenta are distributed is then given by 

d<^n{P;Pl,- ■ ■ ,Pn) —-^gn{pi,P2,- ■ ■ ,Pn)An{pi,P2,- ■ ■ ,Pn) , (30) 

with P = {y/s, 0, 0, 0), where An{pi,p2, ■ ■ ■ ,Pn) is defined in (|^), and where 

, . _ ^P2n-ii^l^^2, ■ ■ ■ ,X2n-4) 

9n[PuP2, . . . ,Pn) - - 3)(logem)2-4 

with 

1 , jPi-l ■ Pi) . 1 1 {Pi ■ Pn) . „ „ 

X2i-?, = log 7- and X2i-4 = ; log -, r , i = 2,... ,n-2 . 

log 4m [Pl-Pn) log 4m [Pi ' Pn) 
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The density obtained before step 4 of Algorithm UA. is invariant under simultaneous scaling 
of all momenta, but the sum of the momenta is not at rest. In order to achieve this, the 
scaling symmetry has to be broken. In order to include the initial-state momenta in the 
density, more symmetries have to be broken: the symmetry under cyclic permutations of 
the momenta, a part of the simultaneous rotation symmetry of the momenta, and finally 
the symmetry under the reflection permutation (1,2,... ,n) {n,n — 1, . . . ,1). 
Let us denote [fc] = k mod n. The cyclic symmetry is broken through 

Algorithm C.2 

1. choose a A; G {0, 2, . . . , n — 1} with relative probability {p^k] ■ Pk+i) 

2. put {pi,P2, • • • ,Pn} {P\l+k-],P\2+k-],- ■ ■ ,P\n+k-]} ■ 

As a result, the function An in the density ( pO]) is replaced by 



n 



{Pn ■ Pi) + {Pl-P2) -\ h {Pn~l ■ Pn 



X (p„ ■pi)An{pi, ... ,Pn 



The factor (p„ ■ pi) in the denominator of An is replaced by the average of all factors. 
Next, part of the rotation symmetry and the symmetry under the reflection permutation 
are broken through 

Algorithm C.3 

1. generate (f G [0,27r] uniformly distributed, and z G [—1, 1 — c] following 1/(1 — 2;) ; 

2. rotate all momenta such that Pi lies along n(2;, (f) . 

If we denote qi := ^y/s{l, 0, 0, 1), this leads to an extra factor 

1 ^ iv^rf 

X 



27rlog(2/c) (qi-pi) 

in the density. The cut-off c can be taken equal to Sq/ (v^p?), where sq should be a cut-off 
on the invariant mass (gi Notice that Pi = {P ■ Pi)/\/s, so that the new factor in 

the density becomes 

1 ^(P-Pi) 



47rlog(2(P -pO/so) (gi-pi) 

In order to include a factor proportional to (p„ -52) in the denominator of the density, 
where ga := |v^(l, 0, 0, -1), Algorithm ^ can be preceded by 

Algorithm C.4 

1. choose with equal probabilities whether to rotate all momenta such that Pi lies along 
h{z, ip), or that Pn lies along —n{z, ip). 
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The total density (|30D becomes then such that An is replaced by 

{P-Pl){Pn-q2) _|_ {qi-Pl){P -Pn) 

^ V log(2(P-pi)/.o) ^ log(2(P.p„)/so) ^ A (n r, r, r, r, n ^ 

47r (p„ ■ Pi) + (Pl ■ P2) H h [Pn-l ■ Pn) 

In the end, we want to obtain all permutations in the momenta {pi,P2, • • • ,Pn, Q'2} of the 
above density. Some of them can be obtained by relabeling, and for the others, which 
change the position of we can do the following. If the cyclic permutation in Algorithm 
|C.2| is chosen with relative probability {p\k'\ ■ Pk+i){Pk+i ■ Pk+i+i) , then the factor {pn-Pi){Pi- 
Pi+i) in the denominator of An is replaced by the average over its cyclic permutations. Then 
we can choose with equal relative probabilities whether to rotate all final-state momenta 
such that Pi lies along n, p^ lies along —n, p^j^i lies along — n, or p„ lies along n. This will 
lead to density that is proportional to 

^n+2(gi, Pl,P2,-- - ,Pi-l,Pi, q2, Pi+l,Pi+2,- ■ ■ ,Pn-l,Pn) ■ 



D Permutations 

For the analyses of the the multi-channel procedure over the antenna densities for the 
different permutations of the momenta, we need an enumeration of the permutations. In 
other words, we need a mapping 

{1,2,... ,n\} ^ Sym^ . 

During a computation, one could, of course, just store an enumeration in the memory 
of the computer, but this costs an amount of memory of 0{n\). Every algorithm that 
delivers all permutations supplies a mapping a priori, but this algorithm has an a priori 
computational complexity of 0{n\). 

We propose the following algorithm as a solution. Any number / G {1,2,... ,n\} can 
uniquely be written in the basis of the lower factorials: I = l+Z2 + ^32!+/43! + - ■ — 1)!, 

where Ik G {0, 1, . . . , — 1}. Let 7^^'^'' denote the z-fold cyclic permutation of the first k 
elements of (1, 2, . . . , n), for example 72 (1, 2, 3, 4) = (2, 3, 1, 4). Then any permutation of 
(1, 2, . . . , n) can be written as 7^2^ ■'7(3^^ ■ ■ ■ 7^"^(1? 2, . . . ,n). This leads to a mapping of the 
required kind, which has a complexity of (9(n^). In a loop delivering all permutations, it 
can easily be reduced to complexity 0{n) at the cost of an amount of memory of 0{n). 
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